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Vortex formation of a Bose-Einstein condensate in a rotating deep optical lattice. 
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We study the dynamics of vortex nucleation and lattice formation in a Bose-Einstein condensate in 
a rotating square optical lattice by numerical simulations of the Gross-Pitaevskii equation. Different 
dynamical regimes of vortex nucleation are found, depending on the depth and period of the optical 
lattice. We make an extensive comparison with the experiments by Williams et al. [Phys. Rev. 
Lett. 104, 050404 (2010)], especially focusing on the issues of the critical rotation frequency for the 
first vortex nucleation and the vortex number as a function of rotation frequency. 
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I. INTRODUCTION 

Ultracold neutral atoms in an optical lattice (OL) are a 
particularly important system for studying a wide range 
of fundamental problems in condensed matter physics [l| . 
When an OL is rotated, the system mimics a lattice sys- 
tem of charged particles subject to a uniform magnetic 
field. This allows the development of versatile quantum 
simulators that can demonstrate various effects caused by 
a magnetic field such as quantum Hall effects 0, 9 • Re- 
cently, two experiments have been reported making use 
of a rotating OL to study quantized vortex dynamics in 
gaseous Bose-Einstein condensates (BECs) [H H| . In this 
system, the vortex pinning parameters can be controlled 
by changing conditions such as the amplitude, lattice con- 
stant, and rotation frequency of the OL. Moreover, the 
realization of a synthetic magnetic field [6j, simulated 
by the Raman process between internal states of atoms, 
opens the possibility of studying a wide range of phe- 
nomena caused by artificial magnetic fields well under 
control. 

Rotating BECs combined with a co-rotating OL ex- 
hibit a rich variety of vortex phases 0412) which have 
two competing length scales: the two lattice spacings, 
one for the OL and the other for the vortex lattice. The 
central role of the OL is to pin vortices at its maxima. 
Tung et al. [H created a rotating square OL using a 
rotating mask, which provided a periodic pinning poten- 
tial that was stationary in the rotating frame associated 
with the presence of vortices. They observed a structural 
crossover from a triangular to a square lattice of vortices 
with increasing potential amplitude of the OL. In a deep 
OL, the condensates are well localized at each potential 
minimum, so that the system enters a regime which can 
be regarded as a Josephson-junction array [13, 0, EH- 
A bosonic Josephson-junction array under rotation (an 
analogue of a Josephson-junction array of superconduc- 
tors under a magnetic field) can have characteristic vor- 
tex patterns with a unit-cell structure which depends on 
a filling factor, i.e., the number of vortices per unit cell of 
the periodic potential. This regime was recently demon- 
strated experimentally by Williams et al. [5j]. 

Motivated by the experiment of Williams et al., we 
study the vortex lattice formation in a rotating BEC 



confined in a deep OL by numerical simulation of the 
Gross-Pitaevskii (GP) equation. We focus on the dy- 
namics of vortices nucleation, which may depend on the 
properties of the OL, finding that the nucleation mech- 
anism is greatly different from the surface instabilit y in 
the case of a harmonically trapped rotating BEC [la - 
lvi\ . An interesting experimental observation is that the 
minimum rotation frequency needed to nucleate a single 
vortex decreases as the amplitude of the OL is increased, 
and eventually falls below the rotation frequency at which 
a single vortex at the center of a condensate is energeti- 
cally stable [f| ■ This implies that vortices were nucleated 
locally at the lattice plaquette. However, our simulation 
reveals that this local nucleation does not occur for the 
experimental parameters used in Ref. We also show 
that inclusion of a Gaussian envelope of the lattice beams 
is important to explain the experimental observation of 
the vortex number dependence on the rotation frequency. 
A preliminary numerical study on vortex nucleation in an 
OL was reported by Yasunaga and Tsubota [l8[ , but the 
case of a deep lattice limit was not discussed and suffi- 
cient understanding of the problem is still lacking. 

The paper is organized as follows. In Sec. HH we de- 
scribe the model of the problem and the parameter set- 
tings for our numerical simulations. We discuss the dy- 
namics of the vortex lattice formation in a BEC under a 
rotating OL in Sec. IIIII and make a detailed comparison 
with experiment in Sec. [IVJ Section [V] is devoted to a 
conclusion. 



II. FORMULATION OF THE PROBLEM 

We study the dynamics of a BEC trapped in an exter- 
nal potential V(r) by employing the dissipative GP equa- 
tion for the condensate wave function \I/(r, t) [iH H9rl2l| : 
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Here, we take a frame rotating with frequency around 
the z-axis. The total external potential V cx t is given by 
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the sum of the harmonic potential Vh and the OL Vol : 
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with r 2 = a; 2 + y 2 \ we shall consider another form of 
Vol later (see Eq. (J7)). . We have denoted that 
L 2 = —ih{xd y — yd x ) is the z-component of the angular 
momentum operator, g — 4irh 2 a s /m is the coupling con- 
stant with s-wave scattering length a s , and 7 is the phe- 
nomenological dissipation parameter. The dissipation 
parameter is introduced to relax the system into the equi- 
librium configuration and is assumed to be 7 — 0.03 in 
the following [ll| [l9[ ; we shall discuss the 7-dependence 
of the results later. 

To rewrite Eq. (fTJ) into a dimensionless form, we use the 
characteristic scales of the harmonic potential as length 
ai 10 = Wh/rruJ±, time , and energy huj±. Since we 
are concerned with the vortex pattern and dynamics in 
the two-dimensional (2D) x-y plane, we reduce Eq. ([T]) 
into the 2D system. We assume that the profile of the 
wave function along the z-axis is approximately uniform 
by considering the system near the z — plane. Then, 
the wave function can be decomposed as &(x, y, z, t) = 
y/N/R z tjj(x, y, t) with the normalization J dxdy\i[i\ 2 = 1, 
where N is the total particle number and R z is a typical 
condensate size along the z-axis, taken as the Thomas- 
Fermi radius. The resulting dimensionless GP equation 
reads 
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with Vol = Vb (sin 2 kx + sin 2 ky) and the dimensionless 
coupling constant m 2 d = 4nNa a /R z . Here, a dimen- 
sionless variable is denoted with a tilde, which shall be 
omitted in the following. 

We use the following parameter values to reproduce the 
experimental setup of Williams et al. [h] : trapping fre- 
quencies (lu±,uj z ) — 27r(20.1, 53.0)Hz, s-wave scattering 
length a s = 5.61 nm, and particle number N = 1.0 x 10 5 
yield ctho = 2.4 fim and R z = 5.6 fim. Then, the cou- 
pling constant becomes w 2 d ~ 1000. The spatial period 
d = 2 /im of the OL used in the experiment corresponds 
to d = n/k = 0.83 in our units. Hence, our free pa- 
rameters are the potential depth Vo in units of hio± and 
the rotation frequency in units of lu±. In the exper- 
iment [5], the range of V is 100 Hz < V < 4000 Hz, 
corresponding to 5 < Vb < 200 in our units. Nucleation 
of vortices was observed by starting from a non-rotating 
condensate loaded into both a harmonic trap and a static 
OL, and then turning on the rotation of the OL. We use 
the ground state solution of Eq. ^ in the presence of 
both the non-rotating harmonic trap and the OL as the 
initial state of the simulations. The numerical scheme to 
solve Eq. (Q} is a Crank-Nicholson method with spatial 
mesh A Xjy = 0.05 and time step A f = 0.0005. Because 
the time development of Eq. Q does not conserve the 



norm of the wave function for 7^0, we treat the chem- 
ical potential /1 as time-dependent and adjust it at each 
time step to ensure normalization, by calculating the cor- 
rection A[i = (At)- 1 \n[f d 2 r\i}(t)\ 2 / f d 2 r\ii(t + At)\ 2 } 



III. DYNAMICS OF A ROTATING BEC IN AN 
OPTICAL LATTICE 

The observation of vortex nucleation and lattice for- 
mation of a harmonically trapped BEC [23j (Vb = 0) 
has been well reproduced by numerical simulations of 
the time-dependent GP equation [l6|, [TtJ ■ The simula- 
tion results clarified the interesting nonlinear dynamics 
where surface wave instability triggers vortex nucleation. 
We first show the overall nonlinear dynamics of vortex 
formation in a rotating BEC subject to an OL to reveal 
the mechanism of vortex nucleation. In this section, a 
rotation with il = 0.6 is suddenly started at t — of the 
simulations. 

Figure [1] shows the time development of the condensate 
density for several values of the potential depth Vq . For 
the shallow OL Vb = 6, the dynamics are similar to those 
without the OL, as shown in Fig. HJa). First, the surface 
of the condensate becomes unstable and generates surface 
ripples that propagate along the surface. Then, the sur- 
face ripples gradually develop into vortex cores, and these 
vortices are pulled in toward the rotation axis to make a 
vortex lattice. The settled vortices are generally pinned 
by peaks of the OL. For V = 20 [Fig. Hp))], although 
vortices are also generated from the surface region, the 
condensate surface is disrupted to form density blobs 
rather than excited to form ripples. The vortices pen- 
etrate inside soon after the turn-on of the rotation. This 
is because quantized vortices accompany the density dips 
(vortex cores) so that they can easily penetrate into the 
condensate through the local density suppression caused 
by maxima of the OL. As the OL becomes much deeper, 
the system enters a regime of the Josephson-j unction ar- 
ray. Here, fractions of the condensate are well localized 
at the potential minima and the overlap of the wave func- 
tion between nearest-neighbor sites is very small. Then, 
the dynamics of the condensate density are completely 
frozen during the overall vortex nucleation process; only 
the phase is a dynamical degree of freedom. A typical 
example for Vb = 100 is shown in Fig. [U(c). Even in 
this deep lattice, the vortices are nucleated from the pe- 
riphery and settle very quickly as soon as the rotation is 
turned on. 

The different behavior in the characteristic vortex nu- 
cleation time can be seen from the development of the 
free energy F = E— ^i — ^(L z ), shown in Fig.QJd), where 
E = J d 2 ri/j*[-{d 2 + d 2 ) + r 2 /2 + V OL + « 2 dMV For 
a shallow lattice, the energy stays constant for a while 
just after the rotation is turned on. During this period, 
the vortices near the surface are prevented from entering. 
After some time, the energy rapidly decreases through 
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(a) j=0 /=5 f = 10 r=20 f = 60 




(b) r = r= 1.6 r = 4 r = 8 f = 20 




FIG. 1: (color online) Time development of the condensate density \ijj\ 2 in the region [—10, 10] x [—10, 10] after the OL suddenly 
begins to rotate at t = with Q = 0.6. The lattice constant of the OL is d = 0.83 (k = 1.2ir) and the height of the OL is Vb = 6 
(a), 20 (b), and 100 (c). In (c), the development of the phase 6 = arg*0 is also shown, where its value changes continuously 
from (dark) to 2n (bright). The vortices correspond to the ends of the branch cuts between the phases and 2n; some of 
them are marked by circles. In (d), time evolutions of the free energy for several values of Vb are plotted as the difference from 
the free energy at t = 0. The physical units are aho = ^/h/muj±, U)J_ , and huj± for length, time, and energy, respectively. 



vortex penetration into the bulk. On the other hand, for 
a deep lattice, this period of constant energy cannot be 
seen. This means that the energy barrier for the vor- 
tex penetration is vanishingly small. Also, the vortices 
equilibrate soon due to the strong pinning effect of the 
OL. We checked that this qualitative feature of the decay 
process was not affected by the values of 7, which sim- 
ply changed the decay time slightly for each numerical 
solution. 

Yasunaga and Tsubota [l8[ found a different kind of 
vortex nucleation dynamic, in which an OL potential 
generates vortex-antivortex pairs inside the bulk region. 
We find that this nucleation dynamic occurs for a rela- 
tively large lattice spacing d. For Vq = 20 but k = O.Stt 
(d = 10/3), for example, a different dynamical feature 
from Fig.QJb) arises as shown in Fig.fSJ some vortices are 
created from the surface as usual, while others arise from 
vortex-antivortex pairs created by peaks of the OL mov- 
ing inside the condensate. The phase profile at t = 0.6 
reveals the pair creation occurring at peaks of the OL [the 



bottom right panel of Fig. [3] , because these vortices have 
circulations and anti-circulations of the phase. Then, the 
anti- vortices quickly disappear by migrating outward and 
combining with other vortices. However, some vortices 
remain in pairs in the condensate and form a vortex lat- 
tice together with vortices coming from the condensate 
surface. We find that similar dynamics occur for different 
values of Vb > 20. This nucleation mechanism is similar 
to that of dragging the superfiow through an obstacle po- 
tential [24j or, more closely, stirring the condensate by a 
circularly moving narrow potential [25j . 

These different dynamical origins of vortex nucleation 
can be understood by noting that the density fluctuation 
is frozen for a deep OL. This condition can be obtained 
when the energy spacing ~ (TrH/d)(2V /m) 1 / 2 of a sin- 
gle well of the OL under the harmonic approximation 
becomes larger than the typical interaction energy ~ ji. 
Then we have d < d c = 27r£^/Vb//i for the density to 
freeze, where £ = h/^/2mJl is the healing length. This 
estimation agrees fairly well with our numerical results 
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FIG. 2: (color online) Time development of the condensate 
density \ip\ 2 after the OL suddenly begins to rotate at t — 
with Q, — 0.6u± for Vo = 20 and k = 0.3-7T (d = 3.3ai 10 ) 
[compare to Fig. Hp))]. The bottom right panel shows the 
phase profile and the inset is an enlarged view of a vortex- 
antivortex pair. 



because d c ~ 0.7 in our parameters (p ~ 83 and £ w 0.1 
for Vo = 100 and Q = 0). The lattice spacing d = 2 /xm 
employed by Williams et al. corresponds to a regime in 
which the local vortex nucleation is not favorable because 
d < d c ~ 3/im for their parameters Vb = 2000 Hz and 
/i = 500 Hz @. 



IV. COMPARISON WITH EXPERIMENTAL 
OBSERVATION 

Here, we compare our results with the experiment by 
Williams et al. [5]. Their significant observations are: 
(i) The minimum rotation frequency of the first vortex 
nucleation decreased with increasing Vo and falls below 
2ir x 1 Hz (Cl = 0.05 in our physical units) above Vb — 
1500 Hz (Vo = 75 in our physical units), (ii) For the deep 
lattice (Vb > /i), the equilibrium vortex number increased 
linearly as a function of the rotation frequency CI. This 
point (ii) indicated that the condensate radius R± was 
not affected by the centrifugal expansion incidental to a 
rotating BEC in a harmonic potential, see Eq. © below. 

In Fig. [3] we show the vortex number N v in the equi- 
librium state as a function of the rotation frequency CI 
for several values of Vb- This result was obtained by the 
dynamical simulation shown in the previous section with 
a sufficiently long time evolution to ensure the equili- 
bration. We can see several features: (i) Naturally, the 
vortex number iV v increases monotonically with increas- 
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FIG. 3: (color online). Equilibrium vortex number as a func- 
tion of Q for several values of the potential depth Vo- We 
count the number of vortices within the radius determined by 
peripheral sites at which the density peak has 10% of the cen- 
tral peak density. The inset shows the behavior of the vortex 
number density taken by the condensate radius at CI = 0, i.e. 
n v = N v /nR ± (n = 0) 2 . 



ing Cl. In contrast to the experiment, it grows faster than 
a linear function for all Vb . (ii) In contrast with Ref. [j| , 
the vortex number N v also increases with increasing Vb 
even at the same rotation frequency. This is because the 
tight confinement by the deep lattice effectively enhances 
the condensate radius R± under the condition of fixed 2D 
particle number [27] . The vortex number thus increases 
with the radius because the lattice spacing of vortices is 
fixed by Cl. However, this property is irrelevant to the 
^-dependence of 7V V , because this 2D artifact only affects 
the radius at Cl = 0. As shown in the inset of Fig. [3J 
when we plot the vortex density taken by the condensate 
radius at Cl = 0, i.e. n v = N v /nR±(Cl = 0) 2 , all the plots 
are almost coincident. This also supports the feature (i), 
namely, the scaling of N v on Cl is independent of Vb- (iii) 
The critical rotation frequency for the first vortex nu- 
cleation is decreased as Vq increases. However, it is still 
above Cl = 0.2 even for the very deep lattice Vb = 100 
in contrast with the experiment [f|. In the following, we 
consider these issues in more detail. 



A. Critical rotation frequency 

Our first concern is the critical rotation frequency of 
the vortex nucleation. To support our results, we calcu- 
late the thermodynamic critical frequency 



Cl r = 



(5) 



that ensures thermodynamic stability of the stationary 
state with a single vortex at the origin [26|. Here, E\ 
and Eq are the total energy for the single vortex state 
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FIG. 4: (color online) Thermodynamic critical frequency fi c 
for the single vortex state as a function of the potential depth 
Vq. The inset shows the corresponding energy difference (nu- 
merator of the right-hand side of Eq. ([5])) and the mean 
angular momentum per atom (denominator of the right-hand 
side of Eq. ©). 



and the non-vortex state, respectively, and (L z ) is the 
mean angular momentum of the single vortex state. Since 
our OL of Eq. ([3]) has a sine form which has a minimum 
at the origin, we replace it with a cosine form with a 
maximum at the origin, which is suitable for calculating 
the energy of the single vortex state pinned at the center. 
Note that this critical frequency is not relevant for the 
actual event of vortex nucleation seen in Fig. [T] because 
the surface instability dominates the vortex nucleation. 
Nevertheless, we can obtain a lower bound of the critical 
rotation frequency because the critical rotation frequency 
fi sur associated with the surface instability is typically 
larger than fi c of Eq. © [TM3]. 

Figure @] shows that the critical rotation frequency fi c 
decreases very slowly with Vo and does not fall below 
0.05 even at Vo = 200. We find that, although the energy 
difference Ei — Eq decreases exponentially with increas- 
ing Vq, (L z ) also decreases together so that their ratio 
51 c does not vanish. For sufficiently large Vq, there is 
little energy cost for a vortex to be positioned at the po- 
tential maximum, because an almost zero density region 
extends between the well separated condensate islands 
at the potential maxima. Correspondingly, the overlap 
of the wave functions between nearest neighbor sites de- 
creases significantly, which leads to the suppression of 
the current flowing between the sites. This reduces the 
atoms contributing to the mean angular momentum. 

Thus, the dynamical simulations discussed in Sec. IIIII 
cannot account for the observation by Williams et al. |5j , 
because the vortices always nucleate from outside of the 
surface and are prevented from invading into the bulk 
for < Q c < f2 sur . Since no vortices were observed in 
a static OL (Q = 0) [Hj], we can exclude the possibil- 
ity of vortex nucleation due to thermal activation like 



the Berezenskii-Kosterlitz-Thouless mechanism [29[ or 
interference between uncorrelated BECs [30(. Once vor- 
tices nucleate, they are expected to have a long life time 
because of the strong pinning effect of the OL, even if 
they are thermodynamically unstable. A more detailed 
study including the effect of the strong fluctuation of the 
OL [5[, or applying a more elaborate model such as the 
projected GP equation (3l| or stochastic GP equation 
[32J to this problem is a challenge for future study. 



B. Rotation frequency vs. vortex number 

The second concern is the linear dependence of the 
vortex number N v as a function of the rotation frequency 
n. For a superfluid rotating in a rigid container with 
radius R, the vortex number N v can be determined by 
Feynman's relation N v = raVlR 2 jh (33[. In the Thomas- 
Fermi limit and if we smooth out the periodicity of the 
density caused by the OL, the radius of a harmonically 
trapped condensate under rotation is given by _Rj_(f2) = 
R±(Q = 0)[1 - (n/cu ± ) 2 }- l/ with v = 3/10 for the 3D 
case and 1/4 for the 2D case [26[. Therefore, the radius 
expands as fl — > uj± and the vortex number diverges as 



mVt 



l 



N v = ^R ± (Q) 2 = —R ± (n = 0) 2 

" (6) 

Contrary to this, Williams et al. observed that for a 
deeper lattice the vortex number increased linearly with 
increasing and did not diverge at fi = cjj_; the vor- 
tex number was nearly coincident with the prediction 
of Eq. © without the factor [1 - (Q / 'u> ±) 2 }~ 2 " , that is 
N v cx fl. This result was expected to be caused by the 
effect of the OL, where the radial expansion may be sup- 
pressed by the small tunneling rate at the peripheral lat- 
tice sites because the number of atoms decreases there. 

We plot the vortex number taken from the numerical 
simulations as well as Eq. ^ in Fig. [SJ To this end, 
we regard the radius R±(£l = 0) as a fitting parameter 
because it is a somewhat arbitrary value for a trapped 
BEC and our interest here is the scaling property of N v 
with respect to Q. For Vo — 20, where the OL plays a mi- 
nor contribution, the vortex number can be well fitted by 
Eq. ([B]), consistent with the experimental observation 
A similar behavior can be seen for Vq = 100 in Fig. OJb). 
Thus, we do not see a linear behavior of the vortex num- 
ber as observed experimentally. 

In order to study the effect of an OL in more detail, we 
need a more realistic model of an OL. Here, we take into 
account the change in Vol due to the Gaussian envelope 
of the lattice beams as 



Vol = -V e 



-2r 2 /t 



[cos 2 (kx) + cos 2 (ky)] (7) 



with beam waist w = 69 /jm [34| . Since the Gaussian 
envelope behaves as 1 — 2r 2 /w 2 + 0(r 4 /w 4 ), it leads to 
an increase in the frequency of the radial harmonic trap 
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FIG. 5: (color online). Comparison of the vortex number 
with Feynman's rule for V = 20 (a) and Vo = 100 (b). The 
numerical results are plotted by filled symbols for the OL 
Eq. © and by empty symbols for the OL with a Gaussian 
envelope Eq. (0. The solid and dashed curves correspond to 
Eq. © and Eq. ©, respectively, with v = 1/4. The fitting 
parameter R±(Q = 0) is 5.0 for both curves in (a) and 6.8 and 
6.1 for the solid and dashed curves in (b). The dotted line in 
(b) represents iV v = R ± (Q = 0) 2 O with R ± (Q = 0) = 6.8. 



as yj uj'j_ + 4Vb I raw 2 . Under this modified potential, al- 
though the dynamical features of vortex nucleation are 
not altered, we can observe a significant reduction of N v 
for large values of Vo, as shown by the empty symbols 
in Fig. [5] This reduction is caused by two contribu- 
tions. The increase in the radial trap frequency makes 
the initial radius R±(tt = 0) smaller than that without 
the Gaussian envelope. We found no visible change of 



R±(n = 0) for V = 20, while R±(Q = 0) becomes about 
10% smaller for Vo = 100. Also, the Gaussian envelope 
suppresses the centrifugal expansion at fl — > u)_\_, which 
allows the condensate radius to remain finite even for 
il > lj±, as observed in the experiment [f|. If we take up 
to the second order of R±/w, Eq. ^ is modified as 



N v = —Rx(Cl = 0) 2 



1 - 



iVo 



4V 



(8) 

which is plotted as the dashed curves in Fig. [5] Although 
only a minor modification occurs for small Vo, it is easy 
to see why the linear-like behavior was observed in the 
range < ft < 1 for large Vq. Our numerical results 
actually appear to follow the line of the linear scaling, as 
seen in Fig. EJb). 



V. CONCLUSION 

We discuss the dynamics of vortex lattice formation 
in a BEC subject to a rotating OL. For a deep OL, we 
found that vortex nucleation occurs in a very different 
way from a usual rotating BEC in a harmonic trap. For 
the lattice spacing d employed by Williams et al. @, 
the motion of the density is frozen and vortices pene- 
trate quickly from the outside through the density chinks. 
Local vortex nucleation, associated with the creation of 
vortex-antivortex pairs in the bulk region, occurs for lat- 
tice spacings larger than the critical value d c . We also 
compare our numerical results with the experimental ob- 
servations of Ref. [5[ , showing that the Gaussian envelope 
of the lattice beam is important in creating a difference 
in the scaling of the equilibrium vortex number on the 
rotation frequency. 
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